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Abstract 


High-order essentially non-oscillatory (ENO) finite-difference schemes are applied to the 
two- and three-dimensional compressible Euler and Navier-Stokes equations. Practical issues, 
such as vectonzation, efficiency of coding, cost comparison with other numerical methods 
and accuracy degeneracy effects, are discussed. Numerical examples are provided which are 
representative of computational problems of current interest in transition and turbulence 
physics. These require both non-oscillatory shock capturing and high resolution for detailed 
structures in the smooth regions and demonstrate the advantage of ENO schemes. 
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1 Introduction 


In the computation of inviscid, compressible flow, the presence of infinitesimally thin shocks 
readily leads to non-linear instability for traditional, unadulterated, linearly stable high- 
order methods. Moreover, regions of strong gradients which have finite thickness but are 
too thin for the grid to resolve may also produce non-linear instability. This is the case, 
for example, for high Reynolds number Navier-Stokes computations in which the shock 
thickness is much smaller than the grid spacing. The standard “cures” are either to add 
explicit artificial viscosity to the numerical method or to employ an upwind-biased scheme 
(which contains implicit artificial viscosity). Such approaches usually have the undesirable 
side effect of loss of resolution, particularly for the small-scale structures in the smooth part 
of the solution. The small-scale features are typically strongly and erroneously damped by the 
artificial viscosity. This problem even afflicts the formally high-order TVD (total-variation- 
diminishing) schemes, since they must degenerate to first order accuracy at smooth critical 
points [13]. Certainly, the most difficult feature to capture is the passage of small-scale 
features through shock waves. 

In many applications, such as typical steady-state aerodynamic CFD, the side effects of 
the artificial viscosity are not particularly worrisome since the main target of the computa- 
tion is the large-scale flow structure and the details of the flow near the shock are not too 
significant. This is decidedly not the case, however, for numerical simulations of transition 
and turbulence. Here the interesting physical phenomena occur on scales much smaller than 
those of the mean flow. As noted by Hussaini and Zang [9], for incompressible flow spectral 
methods have been preferred for these applications since they have the best fidelity for the 
small-scale flow features. However, due to their extreme sensitivity to non-linear instability 
spectral methods have yet to be used for serious investigations of transition and turbulence 
in compressible flow with regions of strong gradients. (They can be used, of course, for 
shock-free flows [4] and for low Reynolds number viscous flows [5] in which thick shocks are 
actually resolved rather than captured .) 

Essentially non-oscillatory (ENO) schemes, first introduced by Harten and Osher [6] and 
Harten, Engquist, Osher and Chakravarthy [7], can achieve uniformly high-order accuracy 
with sharp, essentially non-oscillatory shock transitions. The key idea is an adaptive sten- 
cil interpolation (based on difference tables) which automatically interpolates in a locally 
smoothest region. This strategy provides a strong inhibition towards differencing across 
discontinuities. In [22, 23], Shu and Osher introduced an efficient implementation of ENO 
schemes, using the same adaptive stencil idea but working directly on fluxes and a special 
class of TVD high-order Runge-Kutta type time discretizations. It bypasses the reconstruc- 
tion and Lax-Wendroff procedures in the original ENO schemes. For multi-dimensions, this 
simplification is significant, because the reconstruction in multi-dimensions becomes quite 
complicated [8]. Numerical examples in [22] and [23], especially the examples of shock in- 
teraction with entropy and vorticity waves, for which a good resolution for the detailed 
structures in the smooth region is as important as a sharp, non-oscillatory shock transition, 
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indicate a good potential for ENO schemes in computing compressible Euler and Navier- 
Stokes equations. 

In this paper, we discuss the coding of ENO schemes in [22, 23] to two dimensional general 
geometry (via transformation) and to three-dimensional Euler and Navier-Stokes equations 
of compressible gas dynamics on Cray X-MP and Y-MP. We address the practical issues such 
as vectorization, efficiency of evaluating Newton interpolations, cost comparison with other 
numerical methods, and accuracy degeneracy effects. We then present numerical examples 
which all require non-oscillatory shock capturing and high resolution for rich structures 
in the smooth regions, including two-dimensional shear flows, two dimensional and three 
dimensional homogeneous turbulence, and two-dimensional shock interaction with entropy 
and vorticity waves. 


2 The Navier-Stokes and Euler equations 


For completeness, we document here the three-dimensional, compressible Navier-Stokes equa- 
tions as well as the various transformations that are required for the ENO method in curvi- 
linear coordinates. In terms of the density p , the velocity u = the pressure p, and 

the internal energy E : the Navier-Stokes equations read 

q t -t- f(q)x + g(q)y + h(q) 2 = r(q)* + s(q) y + t(q), (1) 

where 

q = (p,pu,pv, pw, J3)\ f(q) =uq + p(0, 1,0,0, u)‘ 

g(q) = uq + p(0, 0, 1, 0, u)‘, h(q) = u>q + p (0, 0, 0, 1, wf (2) 

and 


- 

? (q) 

= 

/I (0, Tji, T 2 i, T 31 , 0"i) , 


s(q) 

= 

p(0,Ti2,T 2 2,732,<72) , 


t(q) 

= 

P (0,7-13,723, 733, < 73 )*, 

: 

with the components of the viscous 

stress tensor given by 
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r 2 2 

= 
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= 
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4 2 2 
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and 


0i = urn + vTn + WT i3 + 


ct 2 = ur 3 i + ur 33 + wr a3 + 


(7 - 1 )Pr 


(c 2 )„ 

(c 2 ) y 


( 7 - 1 )Pr 

CT 3 = «T 31 + UT32 + WT33 + 7 7 T=- (c 2 ),, 

(7 - l)Pr 


( 5 ) 


Also, \i is the viscosity, 7 is the ratio of specific heats, Pr is the Prandtl number, and 

V = (7-1 )[E - ^p(u 2 + v 2 + w% 


c 2 = — 


if = 


P 

E + p 


( 6 ) 


For implementing ENO schemes with characteristic decompositions, we need the expressions 
of the eigenvalues and the right and left eigenvectors of the Jacobians — , — , The 

di 

eigenvalues for — - are u — c ,11,11,11,11 + c. Its right eigenvectors are the columns of 
dq 


R = 


and the left eigenvectors are the rows of 
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(fe + *) 

(biu + l) 

—biv —b\W 

bi \ 
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0 

2 0 

0 

— 2w 

0 

0 2 

0 

2(1 -b 2 ) 

2biU 

2 i>iu 2b\w 

— 2 bi 


-(hu-l) 

— biv —b\w 

bi / 

bi = 

b 2 = 

7-1 
c 2 ’ 

+ w 2 )bi. 



( 8 ) 


( 9 ) 


dg dh. 

The corresponding expressions for 7^7; and 7^77 are apparent. 
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The Euler equations can be obtained from the Navier-Stokes equations by setting \i = 0. 
The equations for two-dimensional problems are equally obvious. 


In the two-dimensional case, the transformation 

X = x((,r/), y = »((, 7 |), 


( 10 ) 


enables us to treat non-uniform grids or mappings into non-rectangular domains. The Navier- 
Stokes equations become 

q c + fe + &, = *« + £„, (11) 


where 


q = J J q> 

f = J-'m + pfrUtmUy), 

J~ l [Vq + p(0,Tj x ,T)y,Vy], 
J- l (t} x T + ri v s), 

A t (0,T ni T 2 i ) a 1 ) t , 

/i(0,r 1 2,T 2 2,<T2) t , 


g = 

A 

r = 

s = 

r = 

— * 

s = 


with 


n i = g(£*«e + VxUr,) - -(£„«* + v v v v ), 

r 2 l = T12 — -|- TjyUr) -f- £ X V£ T 

r 22 = 7 t (£ v V £ Vy V T)) — n(£x u ( + r ]x U Ti)i 


G 1 = UT n + VT\2 + 


1 


-Mc 2 )t + Vx(c 2 u 


and 


( 7 -l)Pr ] 

<r 2 = UT21 + VT22 + ( 7 _ i)p r [^y( c2 )< + Vv(c 2 ) v ], 


J — £xVy £yVxi 

U = ( X U + £yV, 

V = TJ X U -I- 7}yV. 


( 12 ) 


(13) 


(14) 


The eigenvalues for — are U - c^Jil + 1$, U,U, U + c ^Jil + the right eigenvectors are 
the columns of 

/I c 1 1_ \ 

U - ( X C ~ly U U + {* C 

_ L . v v + & 

\ H - 6 ( c -£ y u + l x v \{u 7 +v 2 ) H + 9 iC ) 


R = 


(15) 
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the left eigenvectors are the rows of 




' (ia + %) ~(M + ^) -(V>_+%) 

2{(y U - L V ) ~ 2 iy 2 ^x 

2(1 - b 2 ) 2 biU 2biv 

{ -(M-%) -(M-4) 

where 6! and b 2 are defined in (2.8); and 

7 £x 


<. = 


JWu}' 

y/G+H' 

(it = Lu + tv v - 

dg 


We again forego the explicit prescription for 


b i 
0 

-26i 

* ) 


(16) 


(17) 


3 Implementing the ENO schemes 


This section should be read in conjunction with [22, 23] for notation, terminology and other 
details of ENO schemes based on fluxes and for TVD time-discretizations of Runge-Kutta- 
type. Here we only summarize several key steps of the algorithms and address practical issues 
such as vectorization, efficiency, cost comparison and the reduction of round-off errors. 

The ENO procedure is applied only to the convection part, i.e., the left-hand-side, of 
Eq (1). The diffusion righthand side of Eq. (1) is approximated by the standard centered 
differences. It is also possible to use ENO-type adaptive stencil interpolation to approximate 
the diffusion terms, but we have not observed any significant differences in our numerical 
tests (typically with small physical viscosity v). Besides simplicity, centered approximations 
also seem more natural for diffusion terms. 


We now summarize the key steps of the algorithm: 

(1) The time marching is implemented by a class of TVD Runge-Kutta type methods [22]. 
For example, the third order case is 

q<’> = q<°> + A(L(q<">), 

q« = ?q (0) + iq W + j A(L(q<‘>), 


4-* ' ’4 

qW = iq<“> + jq (! > + j Aif (q< ! »), 

q <0) = q", 

q n+l = <f>, 


(18) 
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where L is the numerical spatial operator approximating the spatial derivatives in Eq. (2). 
This class of Runge-Kutta methods is labeled TVD because it has been shown [22] that it 
does not increase the total variation of the spatial part under a suitable CFL restriction. 
Also notice that for the third order case, Eq. (18), only three storage levels (two for q, one 
for L) are needed, since can overwrite q^ and q^ can again overwrite q^. 

(2) We thus only need to consider the spatial operator 

L(q) = -f(q)x - g(q)„ - h(q) z + f(q) x + s(q)„ + t(q) z . (19) 


The last three terms are approximated by standard second-order or fourth-order centered 
differences. We use an ENO scheme based on fluxes; hence the first three terms can be 
approximated dimension-by-dimension: for example, when approximating — f(q) x , y and z 
are fixed. The core of the algorithm is then a one- dimensional ENO approximation to, e.g. 
-f(q)x- 

(3) Since f(q) is a vector, we can approximate — f^q)* either component by component, 
or in (local) characteristic directions. In the former case, to obtain non-linear stability by 
upwinding, we write 

f(q) = f t (q) + r(q). (20) 

with 

-w- 1 - 

( 21 ) 


?*(q) = 5(f(q)±<»q), 


df 

where a = max(|u| + c) is the largest eigenvalue in absolute value of the Jacobian — 

- a q 

along the relevant x-line. The decomposition in Eq. (20) guarantees that has posi- 

</<} 

tive/negative eigenvalues; hence upwinding (to be discussed later in this section) is the same 
for all components. For characteristic decompositions, we take Aj+ 1/2 to be some average 
Jacobian at x J+1 / 2 , e.g. the arithmetic mean 


, _ df 

Aj+ * dq 


(22) 


q=3(q,+q i+ i) 


or 


A.-, i = 

J+j 


df 

dq 


-Roe ’ 


(23) 


.(Roe) 


is the Roe average of q ; and q j+1 [17]. We then use the left eigenvectors 


where q v 1 

. 

Rj+ 1/2 in Eq. (8) or Eq. (16) of 1 / 2 to project all relevant quantities (differences around 

Xj+ 1 / 2 ) to the local characteristic fields. A scalar ENO algorithm can then be applied, and 
the result projected back to component space by iZj +1 / 2 in Eq. (7) or Eq. (15). 
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(4) We finally describe the implementation of the scalar ENO approximation of —f(q) x - 
It is written as a conserved flux difference 

-/(?)*!*=*,• w j+i _ ( 24 ) 

where the numerical flux / J+ 1 / 2 approximates /i(xj + i/ 2 ) to a high order with fl(x) defined by 


no#- ( 25 ) 

It is pointed out in [23] that we do not need to construct h(x) explicitly: we simply use the 
difference tables of /(u(x)). If the (undivided) differences of f(u(x)) are defined by 


f\j> 0] = fM, 

= /[? + 1, A: — 1] — f\j, k — 1], Jb = l,-..,r (26) 

where r is the spatial order, then 


r 

/j+i = c(i — j,m) 

m=0 


(27) 


with i being the left-most point in the stencil used to approximate fj+i/ 2 , and c(s,m) being 
defined by 

i J+m *+m 

c(s,m) = (toTTM ^ n (-?)• (28) 

p^l 

The small, constant matrix c is computed once and stored. 

The adaptive stencil determined by the choice of i, the left-most point in the stencil used 
to approximate f j + 1 / 2 . We start with i = j or i = j + 1 according to the (local) wind 
direction (upwinding), and then apply the following 


if (abs(f [i,k]) .gt.abs(f [i-l,k])) i=i-l 


(29) 


for k = 1,- • • ,r. 

Remark 3.1 The code is written in such a way that all the major loops are vectorized by 
default of Cray Fortran. To vectorize Eqs. (27 )-(29) we can either repeat (29) r times (for 
fixed r only) or introduce a temporary one dimensional storage for i to put the short loop (29) 
outside the long loop (27)-(29) for j. To vectorize the characteristic decompositions we have 
to introduce one dimensional temporary storage for the local projection on characteristic 
fields at each Xj. Since we only vectorize the innermost one-dimensional loop, we need just 
17 three dimensional units (10 for two components of q, 5 for L, and 2 work units) for three 
dimensional problems using third order schemes. 
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Remark 3.2 For our current implementation, the componentwise ENO scheme takes 
around 4.5 times as much CPU time as a centered finite-difference scheme with the same 
order of accuracy. A factor of two is due to the flux splitting, Eq. (20). Instead of just 
computing f(q) x we are computing f (q) s + f (q) x ; hence the work is doubled. This is 
the price to pay for implementing upwinding techniques to achieve stability. Another factor 
of two is due to the adaptive stencil procedure, Eq. (29): when these “if” statements are 
removed, the code runs twice as fast. It seems odd that these if statements account for so 
much CPU time since they are all vectorized. The main reason is that since the pattern is not 
uniform from point to point, gathering and scattering are activated by Cray Fortran. These 
procedures are very slow on the Cray. A similar slow-down also exists for TVD schemes. 
ENO schemes using characteristic decompositions take more CPU time: ENO-LF and ENO- 
Roe with entropy fix (see [23] for definitions) take about 3 and 1.7 times, respectively, as 
much CPU time as componentwise ENO schemes. Notice that ENO-Roe is faster than 
ENO-LF because it does not use a flux splitting. See [23] for more details. 

Remark 3.3 We use undivided differences, Eq. (26), and prestored local matrix c, Eq. 
(28), to reduce cost and to reduce the effect of round-off errors. 

4 Transition in a Free Shear Layer 

The numerical examples were chosen to illustrate the ENO method for problems in transition 
and turbulence. We consider flows with gradients which are readily resolved by the grid - 
shocks are either absent altogether or else sufficiently broad (due to low Reynolds number) 
that the usual spectral method gives stable results. Comparisons of the intrinsic resolution 
can therefore be made between the spectral and ENO methods. We then take strong shock 
cases to illustrate the advantage of non- oscillatory high-order methods. 

Unless otherwise indicated, third-order ENO with the third-order Runge-Kutta time dis- 
cretization (18) and fourth-order centered differences for the viscous terms are used. Notice 
that the third-order ENO [23] is actually fourth order in smooth, monotone regions; hence 
for problems with isolated critical points it is fourth-order in L\ norm sense. We use, as 
in [23], the notations ENO-LF (Lax-Friedrichs), ENO-Roe and ENO-Com (componentwise). 

There has been considerable recent interest in the physics of the compressible free shear 
layer and numerical simulations have furnished several interesting results. The numerical 
methods employed have typically been second-order TVD [24, 12], fourth-order MacCor- 
mack [25, 1 6] ,_ or fourth- or even sixth-order compact [10, 19, 2]. Atkins [1] and Sandham 
and Yee [20] have made detailed studies of the performance of TVD schemes on this prob- 
lem. The latter also made comparisions with second-order MacCormack results. Carpenter, 
et al. [2] have compared third-order upwind, fourth-order MacCormack and fourth-order 
compact methods. 

For the particular 2-D examples studied in the present paper, the mean flow is given by 


8 


1 1 Mu I* u mmm i mM a m tunnumumu „ i , H 111, ,1,11 III llllll Mill I IIVIIIIlMlimiilUlilllllii II 1<I 4IIII 1 1 III 1111111111111 


1 

tmT 


where M 0 0 is the ratio 


the hyperbolic tangent profile u 0 = tanh(y), Vq = 0, and po 
in the limit y — > ±oo, of the free stream velocity to the sound speed, and periodic boundary 
conditions are enforced in the streamwise (z) direction. The velocity is non-dimensionalized 

by the freestream velocity, u <*,, lengths by half the vorticity thickness — 2 u 0 n 0 


dy 


the 


density by the freestream density, p OQ} the temperature by the freestream temperature, Too, 
and the pressure by Poo^. The Reynolds number Re = TZo^poo/^oo- The viscosity fi is 
prescribed through Sutherland’s law with a reference temperature of 520°K and the Prandtl 
number is taken to be 0.7. This example is for the temporally evolving free shear layer. A 
forcing term is added to the Navier- Stokes equations in order to make the assumed mean 
flow a steady solution. The computational domain is (0,27r/a) x ( — 00,00), where a is a 
specified wavelength. 


For this problem we present comparisons of ENO with both explicit and compact centered 
difference schemes and with spectral methods. The explicit central difference methods use 
3 points for a second-order approximation and 5 points for fourth-order. The compact 
difference scheme uses a Pade approximation with 5 explicit points and 3 implicit points 
(see [10]). It is formally sixth-order accurate at all interior points, and at points adjacent to 
the boundary, but reduces to fourth-order accuracy at the boundary itself. (The sixth-order 
compact scheme used by Lele [10] reduces to fourth-order accuracy at points adjacent to 
the boundary and third-order accuracy at the boundary itself.) For the spectral calculation 
a Fourier expansion is applied in x and the Cain transformation [3] 


y = —L cot (7?) 


(30) 


is used in the y direction; L is a stretching parameter which is taken between 4 and 10. 
This permits the use of cosine (for p, u , and e) and sine (for u) expansions of the dependent 
variables in the rj direction. 


4.1 Linear Instability 

For the smooth problem we consider first the evolution of a small perturbation, with stream- 
wise wavenumber a — 0.4, from the mean flow at — 0.5 with a Reynolds number of 
100, and the usual 7 = 1.4. The shape of the perturbation is given by the fas test- growing 
eigenfunction of the linearized Navier-Stokes equations at the specified wavelength. The 
amplitude of the perturbation was chosen so that its transverse velocity at y = 0 was 0.1% 
of the freestream velocity. The growth rate for this particular case was 0.127454941. (This 
eigenvalue problem was solved by a spectral linear stability code [11].) For such small am- 
plitudes the non-linear code should produce linear growth of the perturbations for small 
times. Table 1 shows the growth rates produced by central difference methods and compact 
methods after 10 time-steps with a time-step of 0.01. (The measured growth rate was taken 
to be <7 — [log 22(f) - log22(0)]/(2f) where 22(f) = ||tt-ti 0 ||£a + IM|£ a .) In these examples a 
highly-resolved discretization (with 128 points) was employed in y and the specified method 
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method 

N x = 8 

N x = 16 

AT. = 32 

2nd-order central 

-4.26 (-2) 

-6.26 (-3) 

-1.24 (-3) 

4th-order central 

-4.92 (-3) 

-1.83 (-4) 

-1.35 (-5) 

4th-order compact 

2.15 (-4) 

1.25 (-5) 

5.00 (-7) 

6th-order compact 

1.12 (-5) 

1.90 (-8) 

-2 (-8) 


Table 1: Linear Growth Rate Errors for Central-Difference and Compact Schemes at t = 0.1 


method 

OO 

II 

N x = 16 

N x = 32 

II 

2nd-order ENO 
2nd-order ENO-2 

-4.30(-2) 

-3.83(-2) 

-6.30(-3) 

-2.01(-3) 

-1.14(-3) 

-8.28(-4) 

-4.90(-4) 

-5.58(-4) 

3rd-order ENO 
3rd-order ENO-2 

-1.62(-2) 

-8.77(-3) 

-1.05(-3) 

-3.38(-5) 

-1.48(-4) 

-2.20(-6) 

-9.85(-5) 

-2.51(-6) 

4th-order ENO 
4th-order ENO-2 

-4.96(-3) 

-4.21(-3) 

-1.88(-4) 

-1.28(-5) 

-1.54(-5) 

-1.39(-5) 

-9.88(-6) 

-1.78(-5) 


Table 2: Linear Growth Rate Errors for ENO Schemes at t = 0.1 

was applied in x with the streamwise resolution as noted in the table. The table thus pro- 
vides the accuracy achieved as a function of the method and the number of grid-points per 
wavelength. (Even for N x = 4, the error from a spectral discretization in x is already smaller 
than lO -7 .) 

Similar results for ENO methods of second-, third-, and fourth-order are given in Table 
2. In these cases the ENO method was also applied in y, again with 128 points used in 
this direction. Again, the intent was to isolate the discretization errors in x. Except for 
the N x = 64 cases, the results are as one would expect: the accuracy increases with the 
number of grid-points and with the order of the scheme. (The unexpected deterioration of 
the convergence rate for the finest grid is addressed below.) A comparison of the fourth- 
order central- difference results from Table 1 with those of the 3rd-order ENO from Table 2 
indicates that the ENO results are slightly less accurate. This is to be expected since one 
anticipates that in most of the flow this ENO stencil reduces to the fourth-order central one, 
and where it is switched to one-sided it will lose one order in accuracy. 

Notice that the improvement in accuracy obtained from going from 32 to 64 streamwise 
grid-points is less than expected. For the central-difference schemes this occurs on the 10 -6 
level, whereas it occurs an order of magnitude or more earlier for the ENO schemes. (This is 
even more apparent in the ENO results at later times, as evidenced by the data in Table 3.) 
This is due to time-differencing and linearization errors for the central difference methods, 
whereas it is caused by the loss of accuracy when the stencil switches for the ENO method. 
The results labelled by “ENO-2” in the tables are computed by using a factor of 2 to multiply 
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method 

oo 

II 

** 

N x = 16 

N x = 32 

CO 

II 

2nd-order ENO 
2nd-order ENO-2 
3rd- order centered 

-4.26(-2) 

-3.86(-2) 

-4.22(-2) 

-6.37(-3) 

-1.43(-3) 

-6.21(-3) 

-1.26(-3) 

-8.33(-4) 

-1.23(-3) 

-6.20(-4) 

-5.72(-4) 

-5.89(-4) 

3rd-order ENO 
3rd-order ENO-2 

| 

4th- order centered 

-1.59(-2) 

-7.31(-3) 

1.04(-3) 

-1.16(-3) 

1.55(-5) 

3.55(-5) 

-2.41(-4) 

-2.52(-5) 

-3.14(-5) 

-1.89(-4) 

-3.49(-5) 

-3.56(-5) 

4th-order ENO 
4th-order ENO-2 
5th-order centered 

-4.95(-3) 

-4.19(-3) 

-4.92(-3) 

-2.53(-4) 

-5.03-5) 

-2.14(-4) 

-9.35(-5) 

-3.12(-5) 

-4.63(-5) 

-7.79(-5) 

-2.92(-5) 

-4.09(-5) 


Table 3: Linear Growth Rate Errors at t = 1.0 

either the first or the second abs term in Eq. (29), depending upon whether i is greater 
than the left-most point in the centered stencil or not. The effect is to bias the scheme 
towards a centered stencil in smooth regions. This modification of ENO is discussed in 
detail in [21], accompanied by numerical tests on smooth and shocked cases, in response to a 
recent discovery by Rogerson and Meiberg [18] about some accuracy degeneracy phenomena 
of ENO schemes. From the table we can see that “ENO-2” is in most cases comparable in 
accuracy with the corresponding centered schemes, while ENO is usually one order lower, as 
expected from the (unnecessary) switching of stencils in smooth regions. 


4.2 Mach 0.5 evolution 

Next, we present a comparison of these methods for a fully nonlinear problem. The previous 
results were just basic calibration tests (for all the methods). The real purpose of numerical 
simulation codes is to explore nonlinear fluid dynamics. The next example, therefore, is 
a simulation of vortex roll-up and pairing at M 0 0 = 0.5. The initial conditions consist of 
the mean flow plus two linear eigenfunctions: the fundamental with wavenumber oc/ = 0.4 
and amplitude e/ = 0.01 and its subharmonic with wavenumber a, = 0.2 and amplitude 
e, = 0.0001. (In this case the computational domain in x is doubled from that of the 
previous example in order to accommodate the subharmonic.) The initial phases (judged by 
the location of the maximum of the normal velocity perturbation at y = 0) were exactly out 
of phase, a choice which ensures that vortex pairing will occur. 

Figure 1 presents the evolution of the vorticity thickness for third-order ENO on grids of 
size 32 2 , 64 2 , and 128 2 and compares these results with those of a 128 2 spectral calculation. 
(An analysis of the spectral coefficients of the latter coefficients, along the lines discussed 
in [26], indicates that the spectral result is accurate to better than 4 significant digits until 
about t = 125, but that thereafter its accuracy deteriorates rapidly as the vortex roll-up 
produces scales, particularly in the streamwise direction, that are too small for the grid.) 
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The ENO result is clearly converging to the proper answer. A similar comparison is presented 
for the sixth-order compact scheme in Figure 2. The convergence here is more impressive 
than for the ENO scheme, but that is to be expected for this smooth flow. Curiously, the 
spectral coefficients for the compact scheme suggests that it is the transverse resolution which 
is most stressed by the roll-up. 


Figure 3 summarizes, for the third-order ENO method, the evolution of the lowest 4 
Fourier harmonics as represented by the quantity 

E k = f°° (\uk\ 2 + \vk\ 2 )W(y)dy (31) 

J — OO 


where 

n r 2ir/a 

9(z, y,t)f^—dx 

is the k th streamwise Fourier coefficient of the variable q and 



1 

e -(lvl-vc) J 


Is I < y * 
l»l > Vc 


(32) 


(33) 


is a weight function used to confine the region of integration to a finite size. (We used 
y c = 50.) Once again, the numerical results are indicative of convergence. On a 32 2 grid 
the ENO results are perceptibly different from the highly resolved results even for the k = 1 
mode. On a 64 2 grid the worst relative results occur for k — 3. This mode happens to be 
the most sensitive of the 4 to nonlinear interactions. At the start of the calculation only the 
k = 1 and k = 2 modes had non-zero amplitudes. The k = 2 mode is initially forced by 
the self-interaction of the k — 2 mode, which is the dominant mode for the first part of the 
calculation. The k = 3 mode is initially forced by the interaction between the k = 1 and 
k — 2 modes and it is here that the largest errors occur. One heartening result is that the 
A: = 1 mode - the subharmonic - is tracked reasonably well. Atkins [1] observed that there 
could be appreciable spurious generation of this mode by a second-order TVD method. 


Similar data are provided in Figure 4 for the compact scheme. The results for these 
low-order modes are already graphically indistinguishable from the 128 2 spectral results on 
a 64 2 grid for the compact scheme. This, too, is understandable since the compact scheme 
was shown in Table 1 to have an accuracy of better than 1 part in 10 4 with 8 points per 
wavelength, and even the mode k = 4 has 8 points per wave. 

We close the Mach 0.5 results with a plot, in Figure 5, of the pressure contours at 
t = 150 for the ENO method on various grids and for the high-resolution spectral method. 
The similarity of all the results is apparent. 


4.3 Mach 0.9 evolution 

The rationale for the ENO method rests primarily on its behavior in the presence of shocks. 
Indeed, typical central-difference and spectral methods have substantial difficulty for this 
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compressible free layer problem at freestream Mach numbers above 0.70. For this example 
we choose M m = 0.9 and initial conditions consisting of the mean flow plus two linear 
eigenfunctions: the fundamental with wavenumber oif — 0.3 and amplitude tj = 0.01 and 
its subharmonic with wavenumber a, = 0.15 and amplitude e, = 0.001. Furthermore, the 
stretching parameter for the ENO method is here chosen to be L = 10 to provide better 
resolution near the shock waves which eventually develop. 

The evolution of the pressure field for this case is depicted in Figure 6. These plots 
are taken from a computation based on the sixth-order compact scheme. (A 128 3 grid was 
used from t = 0 to t = 75, a 256 grid from t = 75 to t = 100, a 512 x 192 grid from 
t = 100 to i = 106.25, a 768 x 192 grid from t = 106.25 to t = 112.5, and a 1024 x 192 
grid from t = 112.5 to t = 137.5. Spectral interpolation was employed for the requisite grid 
refinements.) By t = 100 vortex pairing has already occurred. The vortex is centered in the 
plotting frame and stagnation points are located on the vertical mid-plane at the streamwise 
edges of the plot. As discussed, for example, by Lele [10], the flow expands away from the 
stagnation points as it goes around the vortex and compresses as it returns to the stagnation 
point. At sufficiently high Mach numbers and for sufficiently strong vortices the compression 
occurs via a shock wave. Shortly after t = 100 in this case a pair of shocks develop — these 
are the so-called “eddy shocklets” [24, 10] - and they grow steadily stronger as the flow 
evolves. These shocks are not the only small-scale feature of the flow, however. As noted 
by Sandham and Yee [20], the flow also develops a very thin region of high strain near the 
stagnation points. 

Although there is a physical viscosity in the flow (in this case the Reynolds number 
Re = 100), the thicknesses of the shock and/or high strain regions may eventually become 
too small for a central difference scheme to handle. Such is the case here even on a 1024 x 192 
mesh. In the presence of unresolved gradients central- difference schemes develop oscillations 
which lead to negative temperatures and an abrupt halt to the calculation. 

Indeed, computations with both the spectral and the sixth-order schemes on 32 , 64 
and 128 2 grids develop severe oscillations and come to a crashing halt between t = 105 
and t = 110. Somewhat curiously, as noted by Sandham (1990, private communication), 
for both methods the oscillations develop first not in the vicinity of the shock but in the 
region of high strain. This is particularly apparent in Figure 7, which shows the evolution of 
the pressure for the spectral method calculation (in which a 128 2 grid was used from t — 0 
to t = 75, a 256 x 128 grid from t = 75 to t = 100, a 384 x 128 grid from t = 100 to 
t = 106.25, a 768 x 144 grid from t = 106.25 to t = 109.375, and a 1024 x 162 grid from 
t = 109.375 to t = 112.0). The computation halted shortly after t — 112 due to negative 
temperatures caused by severe oscillations. Figure 8 is a blow-up of the central region at 
t = 112. Note that there are no oscillations apparent in the vicinity of the shocks. Note also 
that the oscillations are predominantly in the x direction. Indeed, an examination of the 
spectral coefficients reveals that the y direction is quite well-resolved. The price of resolving 
these regions with a non-dissipative central difference scheme can easily be excessive, as the 
present case indicates. 
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The vorticity thickness of the ENO and sixth-order compact methods are provided in 
Figures 9 and 10, respectively. The results for ENO method demonstrate convergence - the 
vorticity thickness on the 128 2 grid is virtually coincident with the compact method result. 
The errors for the ENO method for this Mach 0.9 case are substantially larger than for 
the Mach 0.5 case (see Figure 1). However, this is due primarily to the difference in the 
stretching parameter. For the Mach 0.9 ENO calculations, a weaker transverse stretching 
was used to afford finer resolution in the vicinity of the shock. 

Results for the lowest 4 Fourier harmonics for the two methods are given in Figures 11 
and 12. The performance of the two methods for this diagnostic mimic that for the vorticity 
thickness. 

All of the results reported thus far have been for the ENO method using the characteristic 
decomposition. As noted at the end of Section 3, componentwise ENO is simpler to program 
and is less expensive. Figures 13 and 14 compare the two versions at t = 125 and t — 150, 
respectively. The componentwise results suffer in two respects. First, the shock is more 
diffuse. In fact, at t = 125 the shock is barely visible. Second, there are appreciable spurious 
oscillations. Their character is quite different from the oscillations which afflict the compact 
and spectral results. They are far less regular but are held in check by the nonlinearly 
stable adaptive stencil. Nevertheless, the small-scale flow features of the componentwise 
ENO results are quite unreliable. One must hesitate to use this method for applications 
in which the small scale features are of particular interest, such as transition to turbulence 
problems. 

A comparison of the characteristicwise ENO results at t — 125 with those of the compact 
scheme (Figure 6) reveals that even on a 128 2 grid the ENO method produces a numerical 
shock thickness which is much larger than the actual thickness for this viscous problem. 
Moreover, the absence of a shock on the 64 2 grid appears due to the delayed flow evolution 
(presumably caused by the inherent viscosity of the ENO method) that is apparent in Figure 
9. 

We conclude the results for this problem with Figure 15, which shows the long time 
evolution of a 64 2 characteristicwise ENO calculation based on the Euler equations (but 
starting with the same initial conditions as the Navier-Stokes calculation above). Even on 
this coarse grid, and without the aid of any physical viscosity, the ENO method exhibits 
solid shock-capturing behavior. The numerical solution shows no sign of nonlinear instability 
and spurious small-scale oscillations are absent. 


5 Isotropic Turbulence 

In [15], Passot and Pouquet simulated two-dimensional, compressible, isotropic flows in the 
turbulence regime using a Fourier spectral collocation method. They identified three basic 
regimes: shock-free, weak shocks, and strong shocks. Subsequently, Erlebacher, et al. [5] 
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developed a theory of compressible turbulence that contained a more refined characterization 
of the different regimes of compressible turbulence and contained a useful parametrization of 
initial conditions that permitted precise predictions of the asymptotic turbulence state. The 
direct simulations performed in these studies were limited to quite low Reynolds numbers, 
particularly in the shock regimes. Gibbs oscillations arose whenever the shocks were too thin 
for the grid to resolve. 

Here we perform simulations of compressible isotropic turbulence using both spectral 
and ENO methods. The boundary conditions are periodic in all directions, the velocity 
is normalized by its initial root-mean-square value, the density and the temperature are 
normalized by their mean values, the pressure as for the free shear layer problem, and 
the viscosity is taken to be constant n = 1/150. We compute a low Mach number case 
where the shocks are weak and the spectral method can resolve the full structure with 256 2 
points. Comparisons between different ENO schemes and between ENO schemes and spectral 
methods are furnished for both large-scale and small-scale flow features 

In Figure 16, we show the density and vorticity contours at t = 1 computed with the 
spectral method using 256 2 points. This can be considered to be a resolved solution. Still, 
the vorticity, which involves derivatives for the numerical solution, shows some oscillations. 
In Figures 17 and 18, we show the density and vorticity contours at t = 1 for the spectral 
method and the third order characteristic ENO-LF, respectively, using 64 2 and 128 2 points. 
We can see that the ENO scheme produces non-oscillatory results but the spectral method 
gives noticeable oscillations. For this example, the component ENO-com produces results 
similar to those of ENO-LF. 

In Figures 19, 20 and 21, we show the time history of the average Mach number, the 
maximum Mach number, and the mininum divergence, for the spectral schemes and the third 
order ENO-LF. We can see the convergence of ENO schmes for the former two but not the 
latter. 

In fact minimum divergence is achieved exactly inside the transition regions: to resolve 
it, one has to resolve the full transition regions. The idea of using high order shock capturing 
methods is to resolve essential features in the smooth part of the flow without fully resolving 
the transition regions or shocks. An important topic of numerical tests is to verify whether 
this is achieved. For this example we indeed achieve this as indicated by Figures 19 and 20. 

We then compute a three-dimensional version of this problem. We present in Figure 
22,23,24 the time history of average Mach number, maximum Mach number and minimum 
divergence of velocity. These have characteristics similar to their two-dimensional coun- 
terparts. The minimum divergence time history strongly suggests the presence of three- 
dimensional shocks. However, the grid resolution is not sufficient to clearly bring them out 
by simple flow visualization. 
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6 Shock Interaction With Entropy Waves 


For shock interaction with weak entropy and vorticity waves, some qualitative pictures have 
already been presented in [23, 14], which are similar to those obtained by shock-fitting 
methods [27]. Here we want to do some quantitative studies of wave amplification factors. 
This is relevant to the issue raised in Example 1: if the shock or the rapid transition region is 
not completely resolved, can we still resolve smooth information passing through the shock, 
such as the amplification factor when waves pass through a shock. For this purpose we take 
an entropy wave with a small amplitude so that the linear effects dominate and a comparison 
with linear amplification factor can be made. For third-order ENO with 150 x 20 grid points 
(about 20 points per wavelength), we can already resolve the amplitude amplification factor 
to within 5%. This compares very well against second order MUSCL type TVD scheme with 
the same number of grid points which can only resolve the amplitude amplification factor 
with an error six to ten times as big. Similar comparisons in the one- dimensional case were 
also made (via graphs) in [23]. We remark that this quantitative comparison is important in 
this case, since a major difference between ENO and TVD schemes is that the latter “clips” 
the critical points. If we only compare contours we will not see such sharp differences. 


The details of this problem are as follows. For a pure shock with Mach number M moving 
to the right, we add an entropy wave 

(34) 


— —coafir 

p - p r e Pr 


where )3 T — ^(x cos a r +ysina r ), to the density field at the right of the shock. Here ct T is 
the angle of the vorticity wave with the shock, k T controls the number of waves, and e r is 
the scaled amplitude. In order to enforce periodic boundary conditions in the y direction, 
we take the computational domain to be [0,1] x [0, ^^ — ]. 

The first phase of this computation is aimed at reducing the transients that arise from 
the discrete ENO approximation to the moving shock wave. We run the scheme until the 
shock moves from x — 0.2 to x = 0.8, then shift the data leftwards so that the shock is again 
located at x = 0.2, and repeat this process six times. Then for each fixed x to the left of 
the shock, we perform a Fourier analysis on the entropy to find the amplitude ej, where (34) 
with the the subscripts “r” replaced by “1” denotes the entropy wave to the left of the shock. 
The resulting amplitude ei is then averaged over an x-interval between the wave front and 
the shock, with a length at least one full wavelength. 


The computed amplification factors jf-, together with the linear prediction results, for 
Mach 3, a r = 30°, e? = 0.02, k r = 15, are listed in Table 4. The 6th-order ENO method 
refers to a scheme which is sixth-order in space, and third-order (Runge-Kutta) in time, 
with a reduced CFL number of 0.2. We can observe from Table 4 that, especially for coarse 
grids, higher-order methods indeed produce much more accurate amplification factors than 
low-order methods. It seems, however, that we cannot reduce the error below a certain 
threshhold around 2%. Again, round-of effects might be playing a role since the amplitude 
of the wave is far smaller than the shock strength. 
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method 

N m = 50 

O 

o 

J— 1 

II 

fe? 

N x = 150 

2nd-order MUSCL 

-86% 

-66% 

-42% 

3rd- order ENO 

-47% 

-14.5% 

-6.86% 

6th-order ENO 

-21% 

-8.63% 

-4.74% 

3rd-order ENO-2 

-47% 

-5.98% 

-1.82% 

6th-order ENO-2 

-8.43% 

-6.96% 

-2.38% 


Table 4: Relative Errors in Amplification Factors 

7 Concluding remarks 

ENO schemes based on fluxes and Runge-Kutta type TVD time discretizations, introduced 
in [22, 23] are implemented on Cray 2 supercomputers. Vectorization is realized for all 
inn er loops. Currently the code runs 4.5 times slower than the classical centered difference 
schemes of the same order: a factor of 2 is due to the upwind flux splitting / = / + + /“, 
another factor of 2 is due to the adaptive stencil process. If characteristic decompositions are 
used, the CPU time is increased by another factor of 1.7 to 3. General geometry is handled 
by transformations. Numerical examples include 2D and 3D homogeneous turbulence, shear 
flows, and shock interaction with vorticity waves. ENO schemes show their advantage when 
the solution contains both strong shocks and detailed structures: with a relatively coarse 
grid, where shocks or rapid transition regions are not fully resolved, quantities like mini- 
mum divergence cannot be resolved, but the numerical result is still stable and large-scale 
quantities such as Mach number and amplification factors can be well resolved. 
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Figure 3. Evolution of the 4 lowest Fourier harmonics for the Mach 0.5 free shear layer 
problem using the 3rd-order ENO scheme. 
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Figure 4. Evolution of the 4 lowest Fourier harmonics for the Mach 0.5 free shear layer 
problem using the 6th-order compact scheme. 
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Figure 8. Blow-up of the central region from the spectral results for the free shear layer 
problem at t = 112. 
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gure 9. Evolution of the vorticity thickness for the Mach 0.9 free shear layer problem using 
e 3rd-order ENO scheme. 
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Figure 11. Evolution of the 4 lowest Fourier harmonics for the Mach 0.9 free shear layer 
problem using the 3rd-order ENO scheme. 
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Figure 12. Evolution of the 4 lowest Fourier harmonics for the Mach 0.9 free shear layer 
problem using the 6th-order compact scheme. 
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Figure 15. Pressure field for the Mach 0.9 free shear layer problem from a 64 2 inviscid ENO-2 
calculation. 


















Figures 16: density (left) and vorticity (right) contours for the spectral scheme with 256 2 
grid points. 
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Figures 17: Density (left) and vorticity (right) contours for the spectral scheme with 64 2 
(top) and 128 2 (bottom) points. 
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Figure 19. Time history of average Mach number. 
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Figure 20: Time history of maximum Mach number. 
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Figure 21: Time history of minimum divergence. 
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Figure 24: Time history of minimum divergence, 3D. 
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